Enfoque estadístico del aprendizaje

Trabajo práctico nº1: regresión lineal

 

Datos a utilizar:

Los datos con los que se trabajará en este TP provienen de la 3° Encuesta Mundial de Salud Escolar (EMSE) provistos por el Ministerio de Salud de la República Argentina. Esta encuesta trata sobre temas de salud y hábitos de las personas en la escuela secundaria que pueden impactar en su salud.

Los datasets que se comparten corresponden a un recorte (muestra) del dataset original, luego del tratamiento de valores atípicos e ingeniería de atributos.

Las variables incluidas son:

Variable	Descripción
record	ID de la observación
edad	Edad en años
genero	género de la persona
nivel_educativo	nivel educativo en que se encuentra la persona
altura	altura en centímetros
peso	peso en kilogramos
frecuencia_hambre_mensual	variable categórica que indica la frecuencia con la que la persona considera que pasó hambre en el último mes porque no había suficiente comida en su hogar
dias_consumo_comida_rapida	cuántos días comió en un restaurante de comida rápida en la última semana
edad_consumo_alcohol	edad en qué la persona comenzó a consumir alcohol
consumo_diario_alcohol	cantidad de tragos que la persona habitualmente toma por día
dias_actividad_fisica_semanal	cantidad de días que la persona realizó una actividad física por un total de al menos 60 minutos en la última semana
consumo_semanal_frutas	cantidad de veces que la persona consumió frutas en la última semana
consumo_semanal_verdura	cantidad de veces que la persona consumió gaseosas (al menos un vaso) en la última semana
consumo_semanal_gaseosas	cantidad de veces que la persona consumió snacks/comida salada en la última semana
consumo_semanal_snacks	cantidad de veces que la persona consumió snacks/comida salada en la última semana
consumo_semanal_comida_grasa	cantidad de veces que la persona consumió comidas altas en grasas en la última semana

Consignas:

El objetivo general del trabajo es poder crear una serie de modelos lineales para explicar y predecir el peso de los estudiantes según la información que proporciona la EMSE.

---

 Se cargan las librerías con las que se va a realizar el trabajo:

library(tidyverse)
library(tidymodels)
library(ggplot2)
library(knitr)
library(GGally)

options(scipen=999)

 

Importo el dataset:

encuesta_salud <- read.csv("encuesta_salud_train.csv")

1) Análisis exploratorio

Primero, observo algunos registros del dataset:

encuesta_salud %>% sample_n(5)

Observo las dimensiones del dataset:

glimpse(encuesta_salud)

Rows: 7,024
Columns: 16
$ record                        <int> 502, 26488, 31473, 14154, 36578, 53730, 11892, 35549, 14928, 11754, 40850, 11038, 32832, 31300, 32497, 6961, 39376, 30541, 18309, …
$ edad                          <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, 16, 14, 15, 17, 15, 14, 15, 17, 17, 16, 14, 12, 15, 17, 15, 16, 16, 16, 17, 15, 17, 16, 16…
$ genero                        <chr> "Femenino", "Masculino", "Masculino", "Masculino", "Masculino", "Masculino", "Femenino", "Femenino", "Masculino", "Femenino", "Mas…
$ nivel_educativo               <chr> "2do año/11vo grado nivel Polimodal o 4to año nivel Secundario", "1er año/10mo grado nivel Polimodal o 3er año nivel Secundario", …
$ altura                        <int> 165, 178, 172, 170, 170, 178, 156, 163, 164, 167, 185, 146, 180, 175, 183, 165, 165, 157, 165, 170, 174, 163, 154, 164, 160, 160, …
$ peso                          <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, 100, 33, 62, 70, 80, 60, 47, 50, 50, 70, 75, 55, 64, 59, 86, 51, 71, 49, 70, 65, 65, 53, 5…
$ frecuencia_hambre_mensual     <chr> "Rara vez", "Rara vez", "Nunca", "Nunca", "Rara vez", "Nunca", "Nunca", "Nunca", "Nunca", "Nunca", "Nunca", "Rara vez", "Nunca", "…
$ dias_consumo_comida_rapida    <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1, 3, 0, 0, 0, 0, 1, 0, 6, 0, 1, 0, 2, 0, 2, 0, 0, 5, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0…
$ edad_consumo_alcohol          <chr> "14 o 15 años", "7 años o menos", "Nunca tomé alcohol más que unos pocos sorbos", "14 o 15 años", "16 o 17 años", "8 o 9 años", "1…
$ consumo_diario_alcohol        <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, 5.0, 0.0, 5.0, 0.0, 0.0, 2.0, 1.0, 0.0, 5.0, 0.0, 1.0, 1.0, 0.5, 0.0, 3.0, 3.0, 1.0, 5.0, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1, 4, 0, 1, 6, 5, 7, 3, 0, 7, 5, 2, 2, 4, 2, 7, 4, 7, 0, 6, 4, 0, 3, 2, 0, 2, 2, 3, 0, 4, 2…
$ consumo_semanal_frutas        <chr> "No comí frutas durante los últimos 7 días", "No comí frutas durante los últimos 7 días", "No comí frutas durante los últimos 7 dí…
$ consumo_semanal_verdura       <chr> "4 a 6 veces durante los últimos 7 días", "4 a 6 veces durante los últimos 7 días", "1 vez al día", "4 o más veces al día", "1 a 3…
$ consumo_semanal_gaseosas      <chr> "1 a 3 veces durante los últimos 7 días", "1 a 3 veces durante los últimos 7 días", "4 a 6 veces durante los últimos 7 días", "1 a…
$ consumo_semanal_snacks        <chr> "1 a 3 veces durante los últimos 7 días", "No comí comida salada o snacks en los últimos 7 días", "4 a 6 veces durante los últimos…
$ consumo_semanal_comida_grasa  <chr> "No comí comida alta en grasa en los últimos 7 días", "4 a 6 veces durante los últimos 7 días", "No comí comida alta en grasa en l…

Aquí se puede ver que el dataset cuenta con 7024 observaciones para las cuales se tienen 16 variables (7 de ellas numéricas y 9 categóricas). La variable peso, que es una de las variables numéricas, es la variable a explicar con el resto.

A continuación, observamos si las variables numéricas - y fundamentalmente el peso - están correlacionadas o no. Para ello, se procede a aperturar por el género de cada una de las personas:

encuesta_salud %>%
  select(where(is.numeric), genero, -record) %>%
  ggpairs(mapping = aes(color = genero), title = "Matriz de correlaciones",
          upper = list(continuous = wrap("cor", size = 5, hjust=0.5)), legend = 25) +
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, vjust=0.5), legend.position = "bottom")

En los gráficos, se observa que existe cierta correlación entre entre el peso y la altura (0,57 para las personas de género masculino, 0,44 para las de género femenino. También hay una correlación moderada entre la edad y el peso.

encuesta_salud %>% 
  group_by(genero) %>% 
  summarise(cor = cor(altura, peso))

Luego, observamos qué categorías existen para la “frecuencia de hambre mensual”:

encuesta_salud %>%
  select(frecuencia_hambre_mensual) %>%
  table() %>%
  prop.table() %>%
  sort(decreasing = TRUE)

.
        Nunca      Rara vez Algunas veces  Casi siempre  Dato perdido       Siempre 
  0.689635535   0.206007973   0.083855353   0.011104784   0.005552392   0.003843964 

Para cada una de ellas, observamos cómo se distribuyen el consumo semanal de verdura y el consumo semanal de comida grasa:

tabla = encuesta_salud %>%
  filter(frecuencia_hambre_mensual != 'Dato perdido') %>%
  filter(consumo_semanal_verdura != 'Dato perdido') %>%
  select(consumo_semanal_verdura, frecuencia_hambre_mensual) %>%
  table()
tabla = as.data.frame(tabla)
colnames(tabla) <- c("consumo_semanal_verdura", "frecuencia_hambre_mensual", "q")

order = c("Nunca", "Rara vez", "Algunas veces", "Casi siempre", "Siempre")
type = c("4 o más veces al día", "3 veces al día", "2 veces al día", "1 vez al día", "4 a 6 veces durante los últimos 7 días", "1 a 3 veces durante los últimos 7 días", "No comí verduras ni hortalizas durante los últimos 7 días")


ggplot(tabla) +
 aes(x = frecuencia_hambre_mensual, 
     fill = factor(str_wrap(consumo_semanal_verdura, 20), levels = str_wrap(type, 20)),
     weight = q) +
 geom_bar(position = "fill") +
 scale_fill_hue(direction = 1) +
 coord_flip() +
 theme_minimal() +
 theme(legend.position = "bottom") + 
 theme(legend.title=element_blank()) + 
 labs(x = "Frecuencia de hambre",  y = "Consumo de verdura (proporción)", title = "Consumo de verduras y cantidad de veces con hambre en el mes") + 
 scale_x_discrete(limits = order)

En estos gráficos, se puede observar que aquellos que suelen pasar hambre (siempre o casi siempre) registran un menor consumo de verduras y hortalizas. En los primeros también hay un alto porcentaje de peresonas que comen verduras y/o hortalizas 4 o más veces por día.

2) Modelo inicial

La primera alternativa para modelar el peso, es la siguiente:

\[ E(peso) = \beta_{0} + \beta_{1} altura+ \beta_{2} edad+ \beta_{3} genero + \beta_{4} diasActividadF isicaSemanal + \beta_{5} consumoDiarioAlcohol \]

Ajusto el modelo lineal:

modelo_simple = lm(peso ~ altura + edad + genero + dias_actividad_fisica_semanal + consumo_diario_alcohol, data=encuesta_salud)

tidy_modelo_simple <- tidy(modelo_simple, conf.int = TRUE)
tidy_modelo_simple

El valor de \(\beta_{0}\) nos indica que el peso esperado para una mujer de 0 cm de altura, 0 años, que no realiza actividad física y que no consume alcohol es de -68.92 kg.

El coeficiente que multiplica a la altura \(\beta_{1}\) es 0.65. Significa que ante cada cm adicional de altura el peso promedio aumenta en 0,65 kg manteniendo todas las otras variables constantes

El coeficiente que multiplica a generoMasculino es 1.26, de lo que se deduce que el peso esperado para una persona recién nacida de género masculino de 0 cm de altura, que no realiza actividad física ni consume alcohol es de -67,65 (-68.92 + 1.2626).

El valor de \(\beta_{4}\) nos muestra que ante cada día de actividad física adicional el peso promedio disminuye en 0,08kg (manteniendo todas las otras variables constantes).

El valor de \(\beta_{5}\) nos indica que ante cada trago diario adicional de alcohol el peso esperado de una persona aumenta en 0,007kg (manteniendo todas las otras variables constantes).

Sin embargo, no todas las variables son significativas. Tras realizar el test de significatividad individual para los días de actividad física semanal y para el consumo diario de alcohol se observan p-valores mayores a 0.05 por lo que no se rechaza la hipótesis nula (\(\beta_{k}=0\)). También se puede ver que los intervalos de confianza para ambos coeficientes contienen el 0, lo que nos da la pauta de que no son variables útiles para explicar el peso de una persona.

Por otro lado, el resto de las variables (altura, edad y género) sí resultan estadísticamente significativas para explicar el peso de una persona.

glance(modelo_simple)

El \(R^{2}\) del modelo, el porcentaje de variabilidad del fenómeno que el modelo logra explicar, es de 0,35. Considerando el p-valor, se rechaza la hipótesis nula del test de significatividad global y podemos concluir que al menos una de las variables regresoras sirve para explicar el peso de una persona.

3) Modelo categóricas

Luego, se probará con un modelo que incopore el consumo semanal de snacks y una interacción entre el género y la edad, en lugar de actividad física y consumo de alcohol:

\[ E(peso) = \beta_{0} + \beta_{1} altura+ \beta_{2} edad+ \beta_{3} genero + \beta_{4} consumoSemanalSnacks + \beta_{5} genero.edad \]

Se harán algunas modificaciones para que el nivel basal del consumo semanal de snacks sea “no comí comida salada o snacks en los últimos 7 días:

encuesta_salud$consumo_semanal_snacks <- relevel(as.factor(encuesta_salud$consumo_semanal_snacks), ref = "No comí comida salada o snacks en los últimos 7 días")

modelo_categoricas = lm(peso ~ altura + edad + genero + consumo_semanal_snacks + genero * edad, data=encuesta_salud)

tidy_modelo_categoricas <- tidy(modelo_categoricas, conf.int = TRUE)
tidy_modelo_categoricas

NA

En este caso \(\beta_{0}\) nos indica que el peso promedio de una mujer de 0 años, 0 cm, que no consume comida salada ni snacks es de -64.2 kg.

Las distintas categorías de consumoSemanalSnacks (por ejemplo, consumo_semanal_snacks1 a 3 veces durante los últimos 7 días) reflejan la diferencia de peso entre una persona que no consume snacks vs una persona que consume snacks con la frecuencia que refleja la categoría. No todas las categorías de consumo semanal de snacks resultan significativas, solo algunas de ellas. Vamo

Por otro lado, el coeficiente edad.generoMasculino en este caso refleja que, para los hombres, un año adicional de edad aumenta el peso promedio en 1,22 + 0,39 kg manteniendo constantes el resto de las variables. La variable tiene un efecto significativo para explicar el peso.

Vamos a ver qué porcentaje de la variabilidad explica el modelo

glance(modelo_categoricas)

Por último, considerando que algunas categorías consumoSemanalSnacks no resultan significativa, se va a realizar un test F para evaluar la significatividad conjunta de la variable en su conjunto para explicar al peso.

tidy(anova(modelo_categoricas))

Considerando el p-valor, se concluye que la variable consumo_semanal_snacks es significativa para explicar el peso de una persona.

A continuación, se va a proponer un nueva definición de las categorías de consumoSemanalSnacks: se van a juntar “consume snacks 1 vez al día”, “consume snacks 2 veces al día”, “consume snacks 3 veces al día” y “consume snacks 4 o más veces al día” en una única categoría “consume snacks todos los días”, y se van a volver a estimar los coeficientes.

encuesta_salud$consumo_semanal_snacks_new <- ifelse(encuesta_salud$consumo_semanal_snacks %in% c("1 vez al día", "2 veces al día", "3 veces al día", "4 o más veces al día"), 
                                                    "Consume snacks todos los días", 
                                                    as.character(encuesta_salud$consumo_semanal_snacks))

encuesta_salud$consumo_semanal_snacks_new <- relevel(as.factor(encuesta_salud$consumo_semanal_snacks_new), ref = "No comí comida salada o snacks en los últimos 7 días")

modelo_categoricas_grouped = lm(peso ~ altura + edad + genero + consumo_semanal_snacks_new + genero * edad, data=encuesta_salud)

tidy_modelo_categoricas_grouped <- tidy(modelo_categoricas_grouped, conf.int = TRUE)
#tidy_modelo_categoricas_grouped

tidy(anova(modelo_categoricas_grouped))

Se observa que todas las variables incluídas en el nuevo modelo son significativas (incluídas todas las categorías de consumo semanal de snacks).

A continuación, se va a consultar la variabilidad explicada por el modelo:

glance(modelo_categoricas_grouped)

Con la nueva agrupación no hay una mejora en el \(R^2\) del modelo.

4) Modelos propios y evaluación

A continuación, se sugieren dos posibles modelos adicionales para intentar explicar el peso de una persona:

A: \[ E(peso) = \beta_{0} + \beta_{1} altura+ \beta_{2} edad+ \beta_{3} genero + \beta_{4} consumoSemanalSnacks + \beta_{5} genero.edad \]

B:

\[ E(peso) = \beta_{0} + \beta_{1} altura+ \beta_{2} edad+ \beta_{3} genero + \beta_{4} consumoSemanalSnacks + \beta_{5} genero.edad \]

Realizar 2 modelos lineales múltiples adicionales y explicar brevemente la lógica detrás de los mismos (se valorará la creación y/o inclusión de variables nuevas).

Evaluar la performance del modelo inicial, el modelo categóricas con las categorías redefinidas de la variable consumoSemanalSnacks y los modelos desarrollados en este punto en el dataset de entrenamiento y evaluación (usar dataset “encuesta_salud_test.csv”). La evaluación de performance consiste en comparar la performance en términos del R cuadrado ajustado, RMSE y MAE sobre el set de entrenamiento y en términos de RMSE y MAE sobre el set de evaluación.

¿Cuál es el mejor modelo para nuestro objetivo de predecir el peso? ¿Por qué?

5) Diagnóstico del modelo

Analizar en profundidad el cumplimiento de los supuestos del modelo lineal para el modelo inicial.

plot(modelo_simple)

6) Modelo robusto

Leer el archivo “encuesta_salud_modelo6.csv”. Este último consiste en el dataset original de train con la incorporación de algunas observaciones adicionales que pueden incluir valores atípicos. En particular, observar la relación entre peso y altura ¿Qué ocurre con estos nuevos datos?

Entrenar el modelo inicial con estos nuevos datos y comentar qué se observa en los coeficientes estimados y las métricas de evaluación (R cuadrado ajustado, RMSE y MAE) respecto al modelo entrenado con el set de entrenamiento original. Entrenar un modelo robusto con la misma especificación que el modelo inicial sobre los nuevos datos. Comparar los coeficientes y su performance (RMSE y MAE) respecto al modelo inicial no robusto entrenado en este punto. ¿Qué puede concluir al respecto? Nota: los registros que se suman en este punto son observaciones ficticias que se generaron a partir de observaciones reales del set de datos original.